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ABSTRACT 

Deep spectroscopy of the planetary nebula (PN) NGC 6153 shows that its heavy ele- 
ment abundances derived from optical recombination lines (ORLs) are ten times higher 
than those derived from collisionally excited lines (CELs), and points to the existence 
of H-deficient inclusions embedded in the diffuse nebula. In this study, we have con- 
structed chemically homogeneous and bi-abundance three-dimensional photoionization 
models, using the Monte Carlo photoionization code MOCASSIN. We attempt to repro- 
duce the multi-waveband spectroscopic and imaging observations of NGC 6153, and 
investigate the nature and origin of the postulated H-deficient inclusions, as well as 
their impacts on the empirical nebular analyses assuming a uniform chemical com- 
position. Our results show that chemically homogeneous models yield small electron 
temperature fluctuations and fail to reproduce the strengths of ORLs from C, N, 
O and Ne ions. In contrast, bi-abundance models incorporating a small amount of 
metal-rich inclusions (~ 1.3 per cent of the total nebular mass) are able to match all 
the observations within the measurement uncertainties. The metal-rich clumps, cooled 
down to a very low temperature (~ 800 K) by ionic infrared fine-structure lines, dom- 
inate the emission of heavy element ORLs, but contribute almost nil to the emission 
of most CELs. We find that the abundances of C, N, O and Ne derived empirically 
from CELs, assuming a uniform chemical composition, are about 30 per cent lower 
than the corresponding average values of the whole nebula, including the contribution 
from the H-dcficicnt inclusions. Ironically, in the presence of H-deficient inclusions, 
the traditional standard analysis of the optical helium recombination lines, assuming 
a chemically homogeneous nebula, overestimates the helium abundance by 40 per cent. 
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1 INTRODUCTION 

A long-standing dichotomy in nebular astrophysics is that 
the heavy element abundances derived from optical recom- 
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bination lines (ORLs) are systematically higher than those 
derived from collisionally excited lines (CELs; see Peimbert 
& Peimbert 2006 and Liu 2006 for recent reviews). Temper- 
ature fluctuations (Peimbert 1967), density inhomogeneities 
(Rubin 1989; Viegas & Clegg 1994) and X-ray irradiation 
of quasineutral chemically homogeneous clumps (Ercolano 
2009) have been proposed to explain this discrepancy. How- 
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ever, extensive observations of planetary nebulae (PNe) have 
shown that they are incapable of accounting for all observa- 
tional features (Liu 2006). In order to explain the CEL/ORL 
abundance dichotomy, Liu et al. (2000, L2000 hereafter) 
presented a bi-abundance model. The model assumed that 
the nebula contains two gas components of different abun- 
dances — the diffuse component of "normal" composition 
that dominates the emission of CELs and another consisting 
of cold H-deficient inclusions posited in the diffuse gas that 
dominates the emission of ORLs. The model explains well 
spectroscopic observations from the ultraviolet (UV) to the 
infrared (IR). However, detailed photoionization models are 
still required to test this scenario and address issues such as: 
'How can the nebular physical conditions and its elemental 
abundances be reliably determined in the exiseence of chem- 
ical inhomogeneities?', and 'What is the nature and origin 
of the H-deficient clumps?' To address these issues, this pa- 
per presents the results of three-dimensional photoionization 
modelling for the PN NGC 6153. 

NGC6153 is an archetypal nebula showing a particu- 
larly large ORL versus CEL abundance discrepancy factor 
(adf) of about ten (L2000). A large number of ORLs from a 
variety of heavy-element ions have been detected in this PN. 
NGC 6153 is thus an ideal object to model and to investigate 
the possible causes of the CEL/ORL abundance discrep- 
ancy problem. Early optical studies by Pottasch et al. (1984, 
1986) suggested that NGC 6153 has a very high metallicity. 
The result was further corroborated by the observations of 
the Short Wavelength Spectrometer on board the Infrared 
Space Observatory (ISO-SWS) in the mid-IR ( Pottasch et 
al. 2003). L2000 presented a comprehensive spectroscopic 
analysis of NGC 6153 from the UV to the far-IR and found 
that the C, N, O, and Ne abundances derived from ORLs are 
about ten times higher than those derived from CELs. Using 
simple empirical models, L2000 showed that the discrepancy 
can be explained by assuming the existence of a H-deficicnt 
component embedded in the diffuse gas. Using optical inte- 
gral field spectroscopic observations of three PNe including 
NGC 6153 with the Very Large Telescope Fibre Large Ar- 
ray Multi Element Spectrograph Argus integral field unit, 
Tsamis et al. (2008) provided further evidence for the ex- 
istence of cold H-deficient inclusions and constrained their 
physical sizes to be smaller than ~ 1,000 astronomical units. 
Pequignot et al. (2002, 2003) constructed one-dimensional 
photoionization models that contain two gas components of 
different abundances. The models reproduced satisfactorily 
most of the observed spectral features of NGC 6153. It is 
however difficult to self-consistently treat the diffuse radi- 
ation fields in one-dimensional models and investigate the 
size and spatial distribution of the postulated H-deficicnt 
clumps. This motivated us to carry out further modelling 
using the Monte Carlo photoionization code MOCASSIN (Er- 
colano et al. 2003a, 2005, 2008) capable of dealing with an 
arbitrary nebular geometry and density and abundance dis- 
tributions. 

The paper is organized as follows. In Section 2, previous 
modelling work and observational data used to constrain our 
models are summarized. In Section 3, we introduce the code 
used for the modelling and present the modelling results. In 
Section 4, we discuss the implications and limitation of our 
results, the properties of the H-deficient inclusions and their 



impacts on nebular studies. A summary is given in Section 
5. 



2 OBSERVATIONS AND PREVIOUS MODELS 

In this section, we summarize observations available to con- 
strain the photoionization models of NGC 6153 and in- 
troduce the earlier one-dimensional bi-abundance work of 
Pequignot et al. (2002, 2003), which serves as the starting 
point of our new three-dimensional models. 

2.1 Spectroscopic observations 

Spectra of NGC 6153 from the UV to the far-IR are available 
from the literature (e.g., L2000, Pottasch et al. 1984, 1986, 
2003). For a given wavelength range, there may exist more 
than one set of observations. Due to varying aperture sizes, 
orientations and measurement uncertainties, different stud- 
ies often obtain different fluxes for a given emission line. For 
the current work, we have adopted line fluxes of L2000, who 
present measurements in the optical obtained with ground- 
based facilities, in the UV obtained with the International 
Ultraviolet Explorer (IUE; 1150-3300A), and in the mid- 
to far IR obtained with the Infrared Astronomical Satellite 
(IRAS;7. 6-22.7 /im) and the Long Wavelength Spectrometer 
on board the ISO (ISO-LWS; 43-197 ^m). L2000 obtained 
two sets of optical spectra. The set of spectra with the slit 
oriented along the nebular minor axis allows one to investi- 
gate the spatial brightness distribution of a given line along 
the slit. The other set, obtained by scanning the slit across 
the whole nebular surface, when combined with the total H/3 
flux [logF(H/3) = -10.86 ergs cm^s" 1 ; Cahn et al. 1992], 
yields integrated fluxes of the whole nebula for all the de- 
tected lines. The integrated line fluxes can be directly com- 
pared with the model predictions. Note that L2000 used 
a logarithmic extinction, c(H/3) = 1.30, deduced from the 
Balmer decrement, and the Galactic reddening law of 
Howarth (1983) to deredden the UV, optical and near-IR 
spectra. They adopted a different value c(H/3) = 1.19, de- 
duced from the observed radio continuum flux density and 
the H/3 flux, to deredden the mid- and far-IR spectra. The 
same dereddening algorithm has been used in the current 
study. Here we have added a few mid-IR CELs falling within 
the spectral window of the ISO-SWS in our modelling. Given 
the fact that the size of the ISO-SWS diaphragm was not 
large enough to cover the entire nebula and the telescope 
pointing was offset by 3'.'49 east and 9'.'0 south of the nebular 
centre, aperture corrections have been applied by compar- 
ing the line fluxes obtained with ISO-SWS and those with 
IRAS. All strong emission lines considered in the previous 
modelling (L2000, Pequignot et al. 2003) were included in 
the current work, c.f. the next section for details. 

In addition to line fluxes from the literature, we have 
supplemented the data with new observations using the 
B&C spectrograph mounted on the European Southern Ob- 
servatory (ESO) 1.52 m telescope. The spectra, obtained 
on 2001 May 12, were taken in order to measure several 
near-IR emission lines (He I A10830, [C i] A9850, and [S ill] 
AA9069,9531), important for model constraints. The He I 
A10830 line is particularly important to constrain the helium 
abundance in a bi-abundance model (c.f. Section 4.1). The 
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spectra covered the wavelength range 7700-13254 A, with 
the slit placed along the nebular minor axis. Two long ex- 
posures, each of 1800 s, and one short exposure of 300 s were 
made. The spectra were reduced using the LONG92 package 
in MIDAS following the standard procedures. The spectra 
were bias-subtracted, flat-fielded, cosmic rays cleaned and 
wavelength-calibrated using exposures of a He-Ar calibra- 
tion lamp. They were then flux-calibrated with the IRAI0 
package using wide-slit observations of HST standard stars. 
After dereddening, the spectra were normalized such that 
Hi /(P12)//(H/3) = 0.01106, as predicted by the recombina- 
tion theory, Case B (Storey & Hummer 1995) for T c = 6,080 
K and N c = 3,500 cm -3 . The T c and N c values here are 
taken from L2000. Fluxes deduced from these observations 
are tabulated in Table 1 (see Section 3). 

In an attempt to detect directly the postulated H- 
deficient knots and to constrain their physical sizes, we 
have also carried out long-slit spectroscopy of NGC 6153 us- 
ing the STIS instrument on board the Hubble Space Tele- 
scope (HST). On 2002 June 3-4, two gratings, G430M and 
G750M, were used to cover the spectral range AA3050-5600 
and AA5450-10100, respectively. A long slit of 52" x Of.' 2 was 
placed along the nebular minor axis through the central star 
to cover the brightest parts of the nebula. The spatial reso- 
lution along the slit was 0'.'05/pixel. The standard pipeline 
procedures in iraf/stsdas were used to reduce the data. 
Unfortunately, these data suffered from poor signal-to-noise 
ratio as well as severe contamination by cosmic rays and 
turned out to be of very limited value for the current inves- 
tigation (see Section 4.4). 



2.2 Imaging observations 

Imaging observations of NGC 6153 have been carried out 
with ground-based telescopes (Pottasch et al. 1986) and the 
HST (L2000). Fig. 1 shows the HST images taken under pro- 
gramme GO-8594 (PI Liu) with the WFPC2 camera in the 
F502N and F656N narrow-band filters on 2000 August 130 
Two 600s and two 500s exposures were made in the F502N 
and F656N filters, respectively. The pixel size was 0'.'0455. 
The images were co-added, cosmic-rays removed and cali- 
brated following the standard recipes (e.g. Holtzmann et al. 
1995). They were then dereddened using an extinction coef- 
ficient c(H/3) = 1.30. The two filters trace respectively emis- 
sion of the [O ill] A5007 and Ha lines, although there could 
be some contamination from the nebular continuum emis- 
sion and the [N n] AA6584,6548 lines in the case of F656N. 
We performed aperture photometry and found that Ha and 
[O ill] A5007 line fluxes deduced from HST images agreed 
within 10 per cent with the corresponding values yielded 



1 MIDAS is developed and distributed by the European Southern 
Observatory. 

2 IRAF is distributed by the National Optical Astronomy Ob- 
servatory, which is operated by the Association of Universities 
for Research in Astronomy (AURA) under cooperative agreement 
with the National Science Foundation. 

3 Based on observations made with the NASA/ESA Hubble Space 
Telescope, obtained at the Space Telescope Science Institute, 
which is operated by AURA, Inc., under NASA contract NAS5- 
26555. 



by the optical scanning spectrum of L2000, suggesting in- 
significant contamination to the //5TF656N image from the 
[N n] lines and the nebular continuum emission. 

The azimuthally averaged radial surface brightness dis- 
tribution of Ha was deduced from the HST F656N image 
and plotted in Fig. 2. The profile is used in our models to 
constrain the nebular density distribution. 



2.3 Previous models 

Previous empirical and one-dimensional photoionization 
models of NGC 6153 provide valuable insight and a use- 
ful guide for setting the initial parameters of the new 
3D photoionization models. L2000 constructed several 
parametrized empirical models to match the observed fluxes 
of selected ORLs and CELs. Recombination contributions to 
CELs were considered. They obtained the following results: 

1) Models of a uniform chemical composition and density 
distribution cannot reproduce the large strengths of ORLs; 

2) Chemically homogeneous models with density variations 
predict too weak ORLs and the hydrogen Balmer disconti- 
nuity and a too strong [Ne ill] fine-structure line compared 
to observations; and 3) Bi-abundance models consisting of 
two components, each with a different temperature, density, 
and chemical composition, can account for most of the ob- 
served patterns of NGC 6153. Two different bi-abundance 
models were presented by L2000. One assumed some low- 
density (N c = 700 cm" 3 ) and low-temperature (T c = 500 K) 
H-deficient material embedded in the "normal" nebular gas 
of N e = 5,500 cm" 3 and T c = 9,500 K. In the other model, 
the hydrogen-depleted material was assumed to be fully ion- 
ized and have a high density (2 x 10 6 cm -3 ) and a moder- 
ate temperature of 4,700 K. In both models, the helium and 
heavy element abundances with respect to hydrogen in the 
H-deficient component were respectively 4 and 100 times 
higher than the corresponding values in the diffuse gas. 

Pequignot et al. (2002, 2003) constructed detailed bi- 
abundance photoionization models of NGC 6153 using the 
one-dimensional code NEBU. Their models were composed 
of sectors extracted from a number of spherically symmet- 
ric models, each made of different input parameters except 
the central ionizing source. In the models, a H-deficient 
component with a small filling factor was mixed with a 
normal-composition shell. They found that models contain- 
ing a small amount of dense, cold H-deficient material [M = 
0.0031M Q , N(B.) = 4410 cm -3 , T c = 1390 K] in pressure 
equilibrium with the surrounding gas of "normal" compo- 
sition [M = O.38M , iV(H) = 1170 cm" 3 , T c = 9040 K] 
can well account for the observed fluxes of most ORLs and 
CELs. Compared to the normal gas, the H-deficient compo- 
nent is enhanced in helium and CNONe elements by factors 
of 6 and 100, respectively. Nevertheless, the average elemen- 
tal abundances of the entire nebula are only slightly affected 
by the presence of the H-deficient component, given its small 
mass. Pequignot et al. (2002) pointed out that helium abun- 
dance relative to hydrogen will be overestimated if compo- 
sition fluctuations are not taken into account and it is also 
essential to consider collisional contribution to excitation of 
He I lines to correctly determine helium abundance. 
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Table 1. Comparison of the model predictions and the observations. The observed intensities have been dereddened and normalized 
such that 7(H/3)= 100. Columns (4)-(12) give the ratios of predicted over observed values (departure ratios). 
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6562.8 
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0.08 
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0.95 
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3888.6 
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0.17 


0.95 


0.78 


0.17 


0.95 


He 1 


4471.5 
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1.02 


1.02 


1.02 


0.68 


0.27 


0.95 


0.68 


0.27 


0.95 


He 1 


5875.7 


18.85 


1.01 


1.02 


1.02 


0.67 


0.35 


1.02 


0.68 


0.34 


1.02 
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1.08 


1.09 
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0.72 


0.40 


1.12 


0.72 


0.39 
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He 1 
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4.35 
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1.37 


1.34 
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0.11 


0.97 
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0.11 


0.97 
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1.32 
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0.23 


1.55 
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1.31 


1.21 


0.10 


1.31 


He II 
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1.03 
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0.99 


0.86 


0.13 


0.99 


0.81 


0.21 


1.02 


He II 


4686.0 


12.85 


1.02 


0.97 


0.98 


0.84 


0.13 


0.97 


0.80 


0.21 


1.01 


BJ/H/3 


3646 


0.60 
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0.8 


0.8 


0.75 


0.26 


1.01 
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0.26 
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C 11 


4267.2 


2.42 


0.16 


0.16 


0.16 


0.10 


0.85 


0.95 


0.10 


0.82 


0.92 


C in 


4650+(Ml) 


0.50 


1.11 


0.76 


0.76 


0.41 


0.27 


0.68 


0.40 


0.39 


0.79 


C in 


4187(M18) 


0.08 


1.45 


1.00 


1.00 


0.54 


1.1 


1.64 


0.52 


1.60 


2.12 


N 11 


4041.3 


0.24 


0.19 


0.21 


0.21 


0.15 


0.89 


1.04 


0.16 


0.84 


1.00 


N 11 


4241.8 


0.22 


0.14 


0.15 


0.15 


0.11 


0.66 


0.77 


0.11 


0.62 


0.73 


N 11 


5679.0+(V3) 


0.94 


0.26 


0.29 


0.28 


0.21 


0.72 


0.93 


0.21 


0.69 


0.90 


N in 


4379.0+(M17) 


0.66 


0.60 


0.45 


0.46 


0.30 


0.94 


1.24 


0.30 


1.35 


1.65 


O 11 


4075.0+(V10) 


4.12 


0.15 


0.15 


0.16 


0.10 


0.64 


0.74 


0.10 


0.67 


0.77 


O 11 


4089.3 


0.55 


0.13 


0.13 


0.14 


0.09 


0.83 


0.92 


0.09 


0.86 


0.95 


O 11 


4649.1 


1.39 


0.23 


0.23 


0.24 


0.16 


0.89 


1.05 


0.16 


0.92 


1.08 


O 11 


4651.0+(V1) 


3.68 


0.21 


0.21 


0.23 


0.15 


0.84 


0.99 


0.15 


0.87 


1.02 


Ne II 


3340.0+(V2) 


1.96 


0.09 


0.09 


0.09 


0.08 


0.56 


0.64 


0.08 


0.61 


0.69 


Ne II 


3700.0+(Vl) 


1.07 


0.10 


0.09 


0.10 


0.08 


0.54 


0.62 


0.08 


0.60 


0.68 


Ne II 


4392.0 


0.15 


0.11 


0.11 


0.11 


0.06 


0.56 


0.65 


0.09 


0.61 


0.70 


Mg II 


4481.0 


0.05 








0.66 


0.26 


0.92 


0.67 


0.27 


0.93 




Collisionally excited lines 


(UV and optical) 












[C 1] 


9850.3 


0.17 


0.06 


0.17 


0.34 


0.25 


0.07 


0.32 


0.26 


0.05 


0.31 


C 111] 


1907.0+ 


46.48 


1.22 


1.28 


1.06 


0.90 


0.07 


0.97 


0.87 


0.11 


0.98 


C iv 


1549.0+ 


<30.0 


1.08 


0.61 


0.61 


0.53 


0.00 


0.53 


0.49 


0.00 


0.49 


[N 1] 


5199.0+ 


0.24 


0.02 


0.05 


1.06 


1.16 


0.06 


1.22 


1.12 


0.04 


1.16 


[N 11] 


5754.6 


0.83 


0.30 


0.61 


0.75 


1.03 


0.21 


1.24 


1.08 


0.20 


1.28 


[N 11] 


6583.5+ 


64.20 


0.31 


0.72 


0.99 


1.08 


0.05 


1.13 


1.11 


0.05 


1.16 


N in] 


1744.0+ 


16.30 


0.32 


0.44 


0.33 


0.42 


0.00 


0.42 


0.40 


0.00 


0.40 


[O 1] 


6300.3 


0.78 


0.00 


0.02 


0.74 


0.64 


0.00 


0.64 


0.69 


0.00 


0.69 


[O 11] 


3726.1 


18.9 


0.40 


0.75 


0.89 


0.86 


0.22 


1.08 


0.90 


0.25 


1.15 


[0 11] 


3728.8 


9.69 


0.47 


0.80 


0.91 


0.94 


0.20 


1.14 


0.97 


0.22 


1.18 


[0 in] 


4363.2 


4.19 


0.90 


1.05 


0.93 


0.92 


0.00 


0.92 


0.88 


0.01 


0.89 


[0 in] 


5006.8+ 1189.90 


1.10 


1.19 


1.16 


0.96 


0.00 


0.96 


0.94 


0.00 


0.94 


[Ne m] 


3869.0 


93.91 


0.96 


1.00 


0.95 


1.01 


0.00 


1.01 


0.99 


0.00 


0.99 


[S 11] 


4068.6 


1.03 


0.13 


0.31 


0.52 


0.56 


0.00 


0.56 


0.59 


0.00 


0.59 


[S 11] 


4076.4 


0.35 


0.13 


0.31 


0.52 


0.55 


0.00 


0.55 


0.58 


0.00 


0.58 


[S 11] 


6716.4 


3.16 


0.22 


0.38 


0.57 


0.65 


0.00 


0.65 


0.65 


0.00 


0.65 


[S 11] 


6730.8 


5.26 


0.19 


0.36 


0.56 


0.61 


0.00 


0.61 


0.62 


0.62 


0.62 


[S in] 


6312.1 


1.29 


0.97 


1.33 


1.16 


1.39 


0.00 


1.39 


1.38 


0.00 


1.38 


[S in] 


9068.9 


39.74 


0.86 


1.10 


1.03 


1.09 


0.00 


1.09 


1.09 


1.09 


1.09 


[S in] 


9531.0 


113.68 


0.74 


0.95 


0.89 


0.95 


0.00 


0.95 


0.95 


0.00 


0.95 


[CI m] 


5517.7 


0.55 


1.07 


1.13 


1.01 


1.08 


0.00 


1.08 


1.07 


0.00 


1.07 


[CI m] 


5537.7 


0.70 


0.88 


1.01 


0.90 


0.95 


0.00 


0.95 


0.94 


0.00 


0.94 


[Ar m] 


5191.8 


0.10 


0.83 


0.84 


0.77 


0.94 


0.00 


0.94 


0.91 


0.00 


0.91 


[Ar hi] 


7135.8 


19.00 


1.00 


0.96 


0.97 


1.01 


0.00 


1.01 


1.00 


0.00 


1.00 


[Ar iv] 


4711.4 


2.54 


1.00 


0.97 


0.98 


0.95 


0.00 


0.95 


0.93 


0.00 


0.93 


[Ar iv] 


4740.2 


2.35 


1.00 


0.99 


0.99 


0.96 


0.00 


0.96 


0.93 


0.00 


0.93 



Collisionally excited lines (IR) 



III] 


57.3 /im 


72.43 


1.14 


1.03 


1.02 


0.85 


0.16 


1.01 


0.86 


0.16 


1 


02 


III] 


51.8 fim 


267.60 


0.91 


0.79 


0.85 


0.62 


0.23 


0.85 


0.63 


0.25 





88 


III] 


88.4 ^m 


74.44 


0.84 


0.60 


0.65 


0.49 


0.13 


0.62 


0.49 


0.14 





63 
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Table 1. -continued 



Line 


A (A) a 


-^obs 


S 


El 


E2 


B n 


B c 


B 


K 


B c 


B' 


[0 iv] 


25.9 jim 


118.10 


1.18 


1.18 


1.29 


1.10 


0.16 


1.26 


1.04 


0.29 


1.33 


[Ne II] 


12.8 /im 


23.14 


0.07 


0.11 


0.15 


0.16 


1.00 


1.16 


0.15 


0.68 


0.83 


[Ne in] 


15.6 Jim 


253.50 


1.00 


0.94 


1.00 


0.88 


0.50 


1.38 


0.88 


0.59 


1.47 


[No in] 


36.0 Jim 


32.61 


0.65 


0.60 


0.64 


0.57 


0.18 


0.75 


0.57 


0.21 


0.78 


[S m] 


18.7 jim 


47.28 


1.03 


1.20 


1.18 


1.18 


0.03 


1.21 


1.18 


0.03 


1.21 


[S m] 


33.5 jim 


32.61 


0.81 


0.78 


0.77 


0.80 


0.02 


0.82 


0.80 


0.01 


0.81 


[S iv] 


10.5 Jim 


362.16 


1.09 


0.94 


0.93 


0.87 


0.01 


0.88 


0.87 


0.01 


0.88 


[Ar in] 


9.0 Jim 


56.34 


0.40 


0.37 


0.39 


0.37 


0.01 


0.38 


0.37 


0.01 


0.38 


[Ar III] 


21.8 Jim 


2.69 


0.53 


0.48 


0.51 


0.49 


0.00 


0.49 


0.49 


0.00 


0.49 



refers to the total intensity of the multiplct. 




Figure 1. HST/WFPC2 images of NGC 6153 obtained with the F656N (left panel) and F502N (right panel) narrow-band filters. North 
is up and east to the left. The units of the color bars are 10 -13 ergs cm -2 s — 1 arcsec - 2 . 



3 MODELS 

The modelling is carried out using MOCASSIN (Version 
2.02.12; Ercolano et al. 2003a, 2005), a fully three- 
dimensional Monte Carlo photoionization code that solves 
the radiative transfer self-consistently in an iterative way 
and is designed to tackle problems involving density varia- 
tions and chemical inhomogeneities. The code is therefore 
suitable for the current investigation aiming to study the 
effects of the complicated geometry and possible chemical 
inhomogeneities of NGC 6153. To reduce degree of freedom, 
the ionizing star is assumed to have a blackbody spectral en- 
ergy distribution in our models. A more realistic representa- 
tion of the ionizing star is left to a later work. The nebula is 
approximated by a cubical Cartesian grid of cells of individu- 
ally pre-defined density and chemical composition. For each 
cell, the radiation fields (including the stellar and the diffuse 
component) and physical conditions (electron temperature, 
density and ionic fractions of all elements considered) are 
calculated iteratively. Based on an initial guess of the tem- 
perature and ionization structures, the code calculates ra- 
diation fields, solves the ionization and thermal equilibrium 
equations and then updates the temperature and ionization 
structures. This process is iterated until physical conditions 
are converged in most cells. Finally, the emission spectrum 
of the model nebula is calculated by integrating over the vol- 



ume of the nebula. The code has been applied to construct 
photoionization models of the PN NGC 3918, the H-deficient 
knots in the "born-again" PN Abell30, the Wolf-Rayet PN 
NGC 1501 (Ercolano et al. 2003b,c, 2004) and PN NGC 6781 
(Schwarz & Monteiro, 2006). 

In our models, we have assumed that NGC 6153 has 
an azimuthally symmetrical structure in the x and y plane 
and a reflection symmetry in Z and placed the central star 
at the centre of a corner cell of a 48 x 48 x 48 cubic grid. 
Modelling in this way can effectively save the computational 
time since only one eighth of the nebula is simulated. The 
HST images of NGC 6153 show that the nebula can be ap- 
proximated by an ellipsoidal shell with enhanced emission 
in the equator. The brightest part of the nebula exhibits 
point symmetry with respect to the central star, probably 
implying the existence of an ellipsoidal ring lying along the 
major axis in the south-east to the north-west direction. A 
faint halo can also be seen surrounding the bright nebula. 
In order to better constrain the input parameters, we have 
constructed a series of models with progressively increas- 
ing complexity. These include a spherical model (model S), 
two ellipsoidal models without and with a torus (models El 
and E2, respectively), and a bi-abundance model (model B). 
The spherical and ellipsoidal models are chemically homo- 
geneous. Table 2 gives the adopted model parameters that 
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Table 2. Model parameters. 



Parameter 


S 


El 


E2 


B n 


B c 


B 


L (e36 ergs s _1 ) 


16 


14 


14 


13 


13 


13 


T cf f (kK) 


90 


92 


92 


92 


92 


92 


M(M Q ) 


0.349 


0.269 


0.267 


0.243 


0.0031 


0.246 


iVtot(H) (c56) 


2.63 


2.03 


2.02 


2.07 


0.0081 


2.08 


He 


0.138 


0.138 


0.138 


0.10 


0.50 


0.102 


C (e-4) 


5.4 


4.9 


4.9 


3.2 


177 


3.88 


N (e-4) 


4.6 


4.6 


4.6 


3.8 


150 


4.37 


O (e-4) 


6.82 


6.82 


7.44 


5.63 


440 


7.33 


Ne (e-4) 


1.87 


1.76 


1.90 


1.76 


177 


2.44 


Mg (e-5) 








3.8 


12.1 


3.83 


Si (e-5) 








3.5 


11.3 


3.53 


S (e-5) 


1.75 


1.75 


1.75 


1.75 


5.16 


1.76 


CI (e-7) 


3.0 


2.56 


2.53 


2.35 


10.1 


2.38 


Ar (e-6) 


3.0 


2.75 


3.0 


2.9 


11.5 


2.93 


Fc (e-6) 








1.5 


110 


1.92 


r x a (Hi) 


2.61 


8.91 


362.6 






473 


r z b (Hi) 


2.63 


3.69 


5.275 






6.1 


r x c (He i) 


0.54 


1.83 


73.72 






95.7 


r z d (He I) 


0.54 


0.76 


1.088 






1.24 


r x c (He II) 


74.5 


89.1 


120.5 






90.1 


r z f (He II) 


77.6 


53.7 


59.93 






42.1 


(Tc) (K) 


8,663 


8,825 


8,583 


9,007 


815 


8,892 


t 2 s 


0.008 


0.006 


0.006 






0.014 



a Optical depth along the x axis at the H I photoionization 
threshold. 

b Optical depth along the z axis at the H I photoionization 
threshold. 

c Optical depth along the x axis at the He I photoionization 
threshold. 

d Optical depth along the z axis at the He I photoionization 
threshold. 

° Optical depth along the x axis at the He II photoionization 
threshold. 

f Optical depth along the z axis at the He II photoionization 
threshold. 

8 As definded by Peimbert 1967 

yield the best match to the observations. In Table 1, we com- 
pare the predicted line fluxes with the observations. Further 
details of our models are presented below. 



3.1 Updates to MOCASSIN 

For the purpose of our study, we have made a number of 
updates to MOCASSIN Version 2.02.12, including incorpora- 
tion of He I line emissivities valid at low temperatures, di- 
electronic recombination of the third-row elements of the 
periodic table, recombination contributions to CELs, and 
photoionization depopulation of the He I meta-stable level 
2 3 S. 



3.1.1 He l line emissivities 

The original code uses He I line emissivities at one electron 
density of 10 2 cm" 3 for T c range 5,000-20,000 K, insufficient 
to treat emission from the cold component (T e < 1, 000 K) 
present in our bi-abundance model. In addition, emissivi- 
ties of some He I lines, such as the A10830, are sensitive to 
electron density. Using the calculations of Benjamin et 
al. (1999) and Smits (1996), we have extended the 



He I line emissivities in the code to cover the den- 
sity range 10 2 — 10 6 cm 3 and the temperature range 
312.5—5,000 K. Note that emissivities of a few weak He I 
-lines (A2946, A3615, A4122, A4439, A5049, and A9466) for 
this low temperature range are unavailable. Collisional ex- 
citation from the He I 2 3 S and 2 1 S meta-stable levels is 
insignificant at low temperatures, and is only considered for 
T c > 5,000 K. Using the atomic data of Porter et al. 
(2005), Porter et al. (2007) have provided the most 
recent He I line emissivities. Aver et al. (2010) have 
compared the results of Porter et al. (2007) and Ben- 
jamin et al. (1999). The differences therein are gen- 
erally small compared to the "departure ratios" in 
Table 1. 



3.1.2 Radiative transfer of He I 

For the He I singlets, CaseB recombination is assumed. 
For the He I triplets, the lowest term, 2 3 S, is meta-stable 
and thus considerably populated, leading to significant self- 
absorption and collisional excitation from this level. The 
He I 2s 3 S - 3p 3 P A3888 line suffers the strongest suppression 
due to self-absorption. A fraction of those photons, depend- 
ing on the optical depth, are absorbed and converted into 
_the A7065 line photons. The He I A10830 line is strongly af- 
fected by collisional excitation. In fact, it is mainly excited 
by the collisional process He I 2s 3 S — 2p 3 P. Thus, an accu- 
rate calculation of the population of level 2s 3 S is essential to 
predict reliable fluxes of the He I AA3888, 7065, 10830 lines. 
The 2s 3 S level is populated by recombinations of He + with 
electrons to He triplet levels, and is depopulated by the 
following processes (Clegg and Harrington 1989): (a) colli- 
sional excitation to He singlet levels; (b) radiative decay to 
the ground state; (c) collisional ionization by thermal elec- 
tron impacts; (d) photoionization by UV photons with A < 
2600 A, such as the resonant Hi Lycv photons; (e) decay 
from the He I 2p 3 P state to the ground state following col- 
lisional excitation and the trapping of A10830 photons in 
the nebula. Among these processes, (e) contributes only a 
few per cent of the destruction of 2s 3 S and thus can be ne- 
glected. Processes (a)-(c) have been considered by Benjamin 
et al. (1999) but not process (d), because its effects depend 
on the surrounding radiation fields that are unknown with- 
out detailed modelling. The original mocassin code does 
not consider the destruction of meta-stable helium by the 
photoionization of trapped resonant H I Lya photons, a pro- 
cess that has been shown to be important for compact and 
optically-thick PNe (Clegg & Harrington 1989). This pro- 
cess is introduced in our updated version. With the He I 
2s 3 S level population determined, the self-absorption of the 
He I A3888 line can be treated using a Monte Carlo ap- 
proach. Following the absorption of a He I A3888 photon 
from the 2s 3 S state to 3p 3 P state, it can either scatter back 
or cascade down through the 3s 3 S and 2p 3 P states with 
a probability about 0.9 and 0.1, respectively. In the latter 
case, three photons (A4.3/im, AA7065, 10830) are emitted. 
The self-absorption problem of He I triplets has been inves- 
tigated by Robbins (1968) by solving an integral equation 
of transfer for the ideal case of a uniform sphere expanding 
with a constant velocity gradient. For this ideal case, the re- 
sults deduced from our Monte Carlo approach agree within 
10 per cent with those of Robbins for r(A3888) < 10. 



3.1.3 Di-electronic recombination of S, CI and Ar 

While no reliable di-electronic recombination rates for the 
third-row elements have been published to date, they are 
estimated to be larger than or at least comparable in mag- 
nitude to their radiative counterparts. Considering the fact 
that the rates seem to follow a certain trend with ion- 
ization stage, Ali et al. (1991) suggested that for a given 
ion X l+ , better estimates of the di-electronic recombina- 
tion rates than zero can be obtained by taking the aver- 
age of the rates for the second-row elements C I+ , N l+ , O l+ . 
Dudziak et al. (2000) adopted coefficients empirically cal- 
ibrated with an unpublished model of the PN NGC7027. 
More recently, Ercolano et al. (2004) obtained upper limits 
to the rate coefficients of S 2+ , Cl 2+ and Ar 2+ by modelling 
the PN NGC 1501. In our modelling, we have assumed that 
the di-electronic recombination rates of S, CI and Ar are 
all proportional to those of O, and determined their actual 
values by optimizing the model fit to observations (see Sec- 
tion 4.1). 

3.1.4 Recombination contributions to CELs 

Radiative cascades following recombination (radiative and 
di-electronic) of heavy element ions can partially contribute 
to the emission of CELs, especially for lines from the meta- 
stable levels of low-ionization ions. This process is included 
when calculating the emissivities of CELs from [C 1], [C ill], 
[N 1], [N 11], [N iv], [O 11], [O in], [No iv] and [Ne v]. Here 
the relevant radiative and di-electronic recombination rates 
are taken from Pequignot et al. (1991) and Nussbaumer and 
Storey (1984), respectively, except [N 11] and [O 11], for which 
the total recombination rates calculated by P. J. Storey (pri- 
vate communication) are used. 

3.1.5 Convergence criteria 

MOCASSIN uses the Monte Carlo method to simulate the 
propagation of photons inside a nebula. Thus, one of the 
key component of the code is its convergence criteria. The 
original convergence indicator used in the code is the neutral 
hydrogen fraction, namely a model is deemed as converged if 
the fraction of neutral hydrogen in most cells varies less than 
5 per cent for two adjacent iterations. However, only using 
the fraction of neutral hydrogen as convergence criteria can- 
not guarantee the convergence for other ions, especially for 
cells in the X^/X 1 - 1-1 ^ transition regions. For example, we 
find that in a "converged" model based on the original cri- 
terion, the integrated He 2+ /He ratio may vary 10 per cent 
after additional iterations. To overcome this, we have also 
expanded the convergence criteria to include integrated in- 
tensities of emission lines from helium and heavy element 
ions, such as the He 11 A4686 and [O iv] 25.9 fim, in our 
modelling. 

3.2 Central star 

The properties of the central star are basic input parame- 
ters in the models. Based on the detection of broad O vi 
A3811 and C iv A5801 features, L2000 classified the central 
star of NGC 6153 as a [WC]-PG 1159 H-deficient star. The 
classification is supported by a study of the UV spectrum 
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Table 3. Parameters of the central star. 





^off 


L 








(K) 


(L©) 


(M ) 


(cm s~ 2 ) 


Stoy method 


76,530 


2,034 


0.562 


5.37 


H 1 Zanstra method 


76,260 


2,014 


0.562 


5.37 


He II Zanstra method 


90,090 


3,228 


0.586 


5.47 



of the central star (Pottasch et al. 2003). Using the black- 
body approximation, we have calculated the effective tem- 
perature (T e ff), luminosity, core mass, and surface gravity 
of the central star using the Zanstra method (Zanstra 1931) 
and the Stoy method (Stoy 1933). The results are listed in 
Table 3. In the calculations, we have assumed the distance 
to NGC 6153 to be 1.5 kpc (Pequignot et al. 2003). Both the 
Zanstra and Stoy methods require a measurement of the stel- 
lar continuum flux of the central star. For this purpose we 
have made use of the iTST/WFPC2 F502N image and find a 
stellar continuum flux of 2.75 x 10~ 14 ergs cm~ 2 s _1 A -1 at 
5 010 A after subtracting the nebular background emission. 
The applicability of the Zanstra method also requires the 
nebula to be optically thick, which is probably not the case 
for the H I Lyman continuum given that NGC 6153 is a high- 
excitation class PN and is probably matter bound. Thus the 
Zanstra temperature deduced from the H/3 flux gives only a 
lower limit, while that deduced from the He 11 A4686 line is 
closer to the actual value since the nebula is likely to be opti- 
cally thick in the He 11 Lyman continuum. The Stoy method 
has the advantage that it does not depend on the nebular 
optical depth. However, the applicability of this method re- 
quires the measurement of the total flux of all cooling lines 
from the UV to the far-IR, not an easily achievable task. 
We find that the Stoy temperature is lower than the He 11 
Zanstra temperature, suggesting that the total flux of the 
nebular cooling emission lines might have been underesti- 
mated. 

3.3 Spherical model 

3. 3. 1 Model parameters 

An initially chemically homogeneous spherical model (re- 
ferred as model S hereafter) is constructed. Its density dis- 
tribution is constrained by the azimuthally averaged ra- 
dial surface brightness distribution of Ha deduced from the 
HST/WFPC2 F656N image. We find that a density distri- 
bution given by 

Ps{r,0,$) =N (sxf{r\a 1 ,b 1 ) + (l-s)xf(r\a 2 ,b2)) (1) 
where 

/(r|a ' &) = 4^>) r ^ 3eXp( ^ /b) (2) 

can reasonably reproduce the observed Ha surface bright- 
ness distribution, as shown in Fig. 2. The parameters are: 
ai = 23, 61 = 0.0118, a 2 = 12, b 2 = 0.037, and s = 0.2. 
Note that the shell is actually composed of two density pro- 
files with peaks at r = (121,2 — 3) x 61,2 and FWHMs 
(on, 2 — 3)/6i,2, as illustrated in Fig. 3. In Eqs (1) and (2), 
r ranges from to 0.5, corresponding to a nebular angular 
radius of 0" to 16". A black body spectral energy distri- 
bution for the ionized source is assumed. The model was 
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5 1 

Radius (arcsec) 

Figure 2. Comparison between the observed (dashed line) and 
model-predicted (solid line) azimuthally averaged radial surface 
brightness distribution of Ha of model S. 




5 1 

Radius (arcsec) 

Figure 3. Radial number density distribution of hydrogen atoms 
N(H) of the spherical model (solid line). It is composed of two 
components with different profiles (dashed lines). 



constructed taking the T e ff and luminosity of the central 
star, the nebular elemental abundances and the total mass 
as free parameters. 



3.3.2 Model results 

The spherical model has a total mass of 0.35 Mq and is fully 
ionized by the central star which has a T e g of 90,000 K and a 
luminosity of 4,000 Lq. The nebula is optically thin in both 
the H I and He I ionizing continua, but is optically thick 
in the He II ionizing continuum. The electron temperature 
throughout the nebula is remarkably uniform. We obtain a 
mean electron temperature of 8,663 K, lower than the value 
derived from the [O in] temperature diagnostic lines, and 
a small mean-square temperature fluctuation parameter t 2 
(as defined by Peimbert 1967) of 0.008. 

Fig. 4 shows the radial temperature and density distri- 
butions and the ionization structures of the spherical model. 
In Fig. 4 (also in Figs. 7 and particularly 13), the sharp tem- 
perature dip in the high- ionization region (r ~ 3 — 4") is due 
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Figure 4. Distributions of T c (panel a); N c (panel c); H, He 
(panel b) and O (panel d) ionization structures of the spherical 
model. In panel b, the yellow, blue, green, red, and black lines 
represent those of He 2 " 1 ", He + , He , H + and H°, respectively. In 
panel d, the yellow, blue, green, red, and black lines represent 
those of Q 4+ , Q 3+ , Q 2+ , 0+ and O , respectively. 



to efficient cooling by the strong [Ne v] fine-structure lines 
at 24.2 /im and 14.3 fim. This feature, which is a consequence 
of the low (but not zero) density assumed close to the star 
in the present model, has negligible effect on the model neb- 
ular emission. The decline of temperature for r between 4" 
and 8" is caused by the enhanced cooling of the [O in] lines, 
as the ionic concentration of 2+ increases with r. 

The predicted over observed line flux ratios (hereafter 
'departure ratios') are given in the 4th column of Table 1. 
Most hydrogen and helium recombination lines are well 
matched by the model. The small remaining discrepancies 
between observed and predicted fluxes are within observa- 
tional uncertainties, except for the He I AA7281, 10830 lines. 
The latter two lines with large departure ratios will be fur- 
ther discussed in Section 4. The model however fails to ex- 
plain the strengths of all heavy element ORLs where the 
model predictions are generally about one order of magni- 
tude lower than the observed values. 

The fluxes of CELs from high-ionization species are rea- 
sonably reproduced by the model. For low-ionization species, 
such as [C i], [N i], [O i], [O n] and [S n], the predicted fluxes 
are significantly lower than observed. This can be ascribed to 
the oversimplified density distribution assumed in the spher- 
ical model. The nebula is clearly ellipsoidal and bipolar, and 
appears to be ionization-bound at least in some directions, 
where low-ionization species can survive. We infer that there 
exists a dense torus and present two ellipsoidal models with- 
out and with a torus in the next subsection. 



3.4 Ellipsoidal models 

3.4-1 Model parameters 

To investigate possible effects of the geometry, we have built 
two ellipsoidal models without (model El) and with (model 
E2) a torus. Based on the HST images, we assume that the 
nebula has a semimajor axis of 16" in the polar (z-) direction 
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Figure 5. Distributions of T c (left panel), N(H) (middle panel) and the surface brightness of Ha (right panel) in the x-z plane for 
model El. The contours are of arbitrary units. 




Figure 6. Distributions of T c (left panel), N(H) (middle panel) and the surface brightness of Ha (right panel) in the x-z plane for 
model E2. The contours are of arbitrary units. 



and a semiminor axis of 12" in the equatorial (x- and y-) 3.4-2 Model results 
plane. The density distribution of the ellipsoid is set by 



I a x>\ ^ s d+c|cos(»)|) 
Pe(r - fl ' $)= l + d|cos(g)| 9 ' (3) 

where p 3 was defined in Eq (1). For the ellipsoidal models, 
the parameters a, b, and s are set to the same values as those 
adopted in model S. We introduce three new parameters c, d, 
and g to characterize the density gradient in the equatorial 
plane and the polar direction. The values of c, d, and g are 
constrained by the HST images as well as the observed fluxes 
of low- and high- excitation CELs. Via model fitting, we 
obtain (c, d, g) = (1, 0.8, 0.5) for model El and (0.8, 0.5, 0.5) 
for model E2. A torus with a width of 6 cells and a height of 
1 cell is added in model E2. Its density is assumed to be 2.5 
times higher than the original values of the replaced cells. 
The density distributions of model El and model E2 in the 
x-z plane are displayed in the middle panels of Figs. 5 and 
6, respectively. 



Figs. 5 and 6 display the distributions of the predicted tem- 
perature, hydrogen number density and the surface bright- 
ness of Ha of models El and E2, respectively. Fig. 7 shows 
the model radial distributions of the electron temperature, 
density and ionization structures. In Fig. 8, we compare 
the predicted line fluxes of model E2 with observations. We 
find that both models El and E2 give a better fit to the 
[N n], [O n], and [S n] lines compared to model S thanks to 
the density enhancement in the equatorial plane. Model E2 
yields a better fit for the [C i] and [N i] lines than model 
El, owing to the presence of a dense torus. For the remain- 
ing lines, models El and E2 give similar results of model 
S. Although the ellipsoidal models show improvements in 
reproducing CELs from low-ionization ions, like model S, 
they fail to reproduce the strengths of heavy element ORLs, 
strongly indicating that chemical inhomogeneities have to 
be invoked. 
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Mod. E1 : X direction 



Mod. E 1 : Z direction 



Mod. E2: X direction 



Mod. E2: Z direction 
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Figure 7. T c and N c distributions and H, He and O ionization structures in the two ellipsoidal models. Panels (a) — (d) and (c) - 
(h) are for the x, z directions in the ellipsoidal model El, respectively; Panels (i) - (1) and (m) - (p) are for the x, z directions in the 
ellipsoidal model E2, respectively. In panels (c), (g), (k), and (o), the yellow, blue, green, red and black lines represent the fractions of 
He 2+ , He+, He , H+ and H°, respectively. In panels (d), (h), (1), and (p), the yellow, blue, green, red, and black lines represent the 
fractions of 4+ , 3+ , 2+ , + and O , respectively. Note that the sharp features near 10'.'5 in panels (i), (j), (k) and (1) are produced 
by the dense torus. 



2.0 



1.5 



Q 
o 



1.0 



0.5 



0.0 



(a) 



CD OJ CD CO 

co r-- co 

_ -tf CO CO 

-r — — <D 

CD ^ 



illnill 



III 



2.0 



1.5 



1.0 



0.5 



0.0 



(b) 



LO CD 
= £0 



cn en _ in eg m -■- 



- J ^CD — 
= CO I" CO W 
qj CO 1X5 o 2 



Figure 8. Comparison of line fluxes predicted by model E2 and observations for ORLs (left panel) and CELs (right panel). 
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3.5 Bi-abundance model 

3. 5. 1 Model parameters 

For the bi-abundance model, we assume that there are some 
metal-rich knots (cells) embedded in the diffuse nebula of 
"normal" abundances. The density distribution of the diffuse 
nebula is inherited from model E2, but with a slightly re- 
duced mass. The spatial distribution of the metal-rich cells is 
set randomly by the probability density function /(r|a3,fe3) 
(see Eq (2)), where r ranges from to 0.3125 (i.e. 0-10"), 
ci3 = 22 and 63 = 0.013, as constrained by the surface bright- 
ness distribution of the O 11 A4649 ORL along the nebular 
minor axis (see Fig. 19 of L2000). In Fig. 9, we compare the 
profile of the predicted O 11 surface brightness with obser- 
vations. In the comparison, we have included the effects of 
seeing and convolved the predicted surface brightness dis- 
tribution with a Gaussian function of FWHM 3'.'2. Posi- 
tions of individual knots are set by Monte Carlo simulations. 
Their three-dimensional spatial distribution is illustrated in 
Fig. 10. Fig. 11 shows the distributions of number density of 
knots and of the hydrogen atoms of the diffuse gas. In order 
to reduce the random uncertainties, we have generated four 
samples of knots. Line fluxes yielded by the four samples 
are in good agreement (within 10 per cent). They are then 
averaged to compare with the observations. The model pa- 
rameters are given in Table 2, where B n and B c represent re- 
spectively components of "normal" composition and of cold 
metal-rich knots, and B is the sum of the two components. 
Note that for the normal component, we have fixed the he- 
lium abundance to 0.1 (as adopted by Pequignot et al. 2003). 
As pointed out by Pequignot et al., observations of helium 
recombination lines in the optical alone are not sufficient to 
decompose the helium abundance in the two components. 
For the H-deficient knots, the element abundances heavier 
than neon are assumed to be the corresponding solar val- 
ues in mass, considering that PNe originate from low- and 
intermediate- mass stars that can not process the third-row 
elements. 



Figure 10. Three-dimensional spatial distribution of metal-rich 
knots (cells) in the bi-abundance model. The units of the axes are 
arcsec. 




Radius (arcsec) 

Figure 11. Radial distribution of the volume number density 
of metal-rich knots given by the analytic (solid line) and Monte 
Carlo simulations (asterisks). The dashed and dash-dotted lines 
represent the density distribution of hydrogen atoms in the diffuse 
nebular gas along the x and x axes, respectively. 

3.5.2 Model results 

Fig. 12 displays the predicted structure of T e and N e . The ra- 
dial distributions of temperature, density and the ionization 
structures in four directions across the four metal-rich knots 
marked in Fig. 12 are shown in Fig. 13. An inspection of 
the two figures shows that the knots have been cooled down 
to very low temperatures (~ 800 K), and create shadows 
that block the ionizing UV radiation field from the central 
star, leaving low ionization regions behind them. As in mod- 
els S, El and E2, there is a temperature dip at r ~ 3 — 5". 
Fig. 14 displays the predicted monochromatic images of H/3, 
[O III] A5007, and O n A4649. The figure clearly shows that 
the [O in] CEL originates from the "normal" component, 
while the O 11 ORL is dominated by knots. In spite of 
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their small number, metal-rich knots have a non-negligible 
contribution to H/3 emission, owing to the enhanced emis- 
sivity at low temperatures and high densities. In Fig. 15, 
we compare the predicted line fluxes with the observations, 
demonstrating that the bi-abundance model has significant 
improvements in reproducing the strengths of heavy element 
ORLs and the hydrogen Balmer discontinuity. In addition, 
the model yields a much better fit to the measured flux of 
the [Ne n] 12.8 /im line, which is underestimated by all the 
chemically homogeneous models. This is due to the fact that 
the metal-rich knots are enriched in Ne and have a large Ne + 
ionic concentration. 

Model B generally supports the results of the 
one-dimensional models of Pequignot et al. (2002, 
2003) in terms of properties of the H-deficient inclu- 
sions such as the electron temperature, density, mass 
and chemical abundances. Those values in model B 
(T e = 1,390 K, N(E) = 4,410 cm" 3 , M = 0.0031M© 
and oxygen enrichment factor of about 100) are 
comparable to those given by the one-dimensional 
models (T c = 815 K, iV(H) = 4,000 cm" 3 , M = 
O.OO31M0 and oxygen enrichment factor of about 
80). Beyond constraining the properties of the H- 
deficient inclusions in a self-consistent manner, the 
three-dimensional model enables to investigate sizes 
and spatial distributions of the H-deficient inclu- 
sions which one-dimensional models can not. De- 
tailed line by line comparison between model B 
and the one-dimensional models by Pequignot et al. 
(2002, 2003) is not easy, due to different geomet- 
ric configurations, density structures, spectral en- 
ergy distributions of the ionizing star and atomic 
data used. We do not claim superiority of pre- 
dicted fluxes of model B over the one-dimensional 
models. To understand how three-dimensional bi- 
abundance models advance the field of modelling of 
emission line gaseous nebulae, simple bi-abundance 
benchmark models are needed. Then the compar- 
ison between three-dimensional model results and 
one-dimensional model results will be much easier 
to carry out and interpret, but this is beyond the 
scope of the current paper. 



4 DISCUSSION OF THE BI-ABUNDANCE 
MODEL 

4.1 Uncertainties and discrepancies 

Uncertainties in our models include errors in the atomic 
data, in the treatment of the radiation transfer and un- 
certainties in the adopted line fluxes. The estimates of the 
di-electronic recombination rates (aa) for S, CI and Ar are 
rather rough. We have assumed that a<i(S, CI, Ar) = rj(S, 
CI, Ar) x ad(0), where the constants ??(S, CI, Ar) are ob- 
tained by fitting the observations. Only two [CI ill] lines 
are detected, insufficient to place a constraint on the 77(d) 
value. We have thus assumed rj(Cl) = 1. S and Ar have a 
number of lines detected from several ionization stages. We 
obtain r?(S, Ar) = (0.76, 6.8), (0.45, 4.0), (0.35, 4.0) and 
(0.35, 4.0) from the models S, El, E2 and B, respectively. 
In the bi-abundance model, the di-electronic recombination 



coefficients of S and Ar are respectively 0.32 and 14 times 
the corresponding radiative recombination values. Our esti- 
mates of a d (S, Ar) differ from those obtained from the mod- 
elling of NGC 7027 (Dudziak et al. 2000) by a factor of three. 
On the other hand, these uncertainties in the di-electronic 
recombination rates of S and Ar are unlikely to have a large 
impact on our model predictions for other lines, considering 
the low abundance of Ar and the fact that [S 11], [S 111] and 
[S iv] lines have similar emissivities, and consequently ther- 
mal structure is hardly affected by the uncertainties in the 
ionization structure of S. 

The bi-abundance model underestimates the flux of 
C in multiplet Ml, but overestimates that of M18. This 
is probably caused by errors in the atomic data. The ex- 
isting calculations of the radiative (Pequignot et al. 1991) 
and di-electronic (Nussbaumer & Storey 1984) recombina- 
tion coefficients for C ill may not be applicable to plasma 
at very low temperatures. Further theoretical work is ur- 
gently needed to study the behavior of recombination lines 
in extremely cold plasma. We note that because of a larger 
di-electronic contribution the emissivity of lines from C ill 
multiplet M18 is much more sensitive to electron tempera- 
ture than that of Ml. Should it not be the case, the model 
would have yielded a much better fit to both Ml and M18. 

The bi-abundance model overestimates the flux of the 
He 1 A10830 line by 21 per cent. The He I A10830 line is 
crucial for the bi-abundance model in that, unlike other 
He 1 lines, it is predominantly collisionally excited, and thus 
has the potential to constrain the helium abundance in the 
hot component of "normal" composition. Its usage is how- 
ever hindered by the complexity in the He I 2s 3 S meta-stable 
level population and the excitation and radiative transfer 
of the He I A10830 line. Observations of several PNe have 
shown that the measured He I line ratio 7(A10830)//(A5876) 
is lower than the theoretical prediction (see, e.g., Peimbert & 
Torres-Peimbert 1987a,b). Possible explanations include the 
destruction of He I A10830 photons by internal dust (King- 
don & Ferland 1995), and the photoionization depopulation 
of the He I meta-stable level by Hi Lya photons (Clegg & 
Harrington 1989). The later process has been considered in 
our model. However its inclusion reduces the predicted flux 
of the He I A 10830 line only slightly. Additional destruction 
mechanisms may be needed to reconcile observations with 
theory. 

Given its high sensitivity to electron temperature, the 
He 1 A7281 line serves as an important diagnostic to probe 
the physical conditions in the cold H-deficient knots (Zhang 
et al. 2005a). Although improved compared to the chemi- 
cally homogeneous models, our bi-abundance model overes- 
timates the He I A7281 line by 55 per cent. The discrepancy 
is probably caused by departure from the CaseB assump- 
tion. Under Case A recombination, the predicted flux of the 
He 1 A7281 line is about a factor of two lower. L2000 and 
Liu et al. (2001) find that the observed fluxes of the 2s 1 S- 
np P and 2p P-ns 1 S series are systematically lower than 
predicted by Case B recombination. Similar discrepancies 
are also found in H 11 regions, such as the Orion Nebula 
(Porter et al. 2007; Blagrave et al. 2007) and 30 Doradus 
(Tsamis & Pequignot 2005). Porter et al. (2007) pointed 
out that although the escaping of the He I Lyman 
photons could cause the discrepancies, it is not con- 
sistent with the fact that the He I A6678 line flux 
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Figure 12. Predicted JV c (left panel) and T c (right panel) distributions in the bi-abundance model B. 
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Figure 13. Distributions of T e , N c and H, He and O ionization structures along a radial direction that passes through the marked four 
knots in Fig. 12. Panels (a) - (d) and (e) — (h) are for the Knots 1 and 2 located in the xy plane, respectively; Panels (i) - (1) and (m) 
— (p) are for the Knots 3 and 4 in the yz plane, respectively. In panels (c), (g), (k) and (o), the yellow, blue, green, red, and black lines 
represent the fractions of He 2+ , He"*", He , H+ and H°, respectively. In panels (d), (h), (1) and (p), the yellow, blue, green, red, and 
black lines represent the fractions of O 4 ^, O a+ , 2+ , + and O , respectively. Note that in panels (b) and (f), the peaks at ~ 8" and 
~ 11" correspond respectively to the knot and the torus. 
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Figure 14. Projected monochromatic images of H/3 (left panel), [O III] A5007 (central panel) and O II A4649 (right panel) pre- 
dicted by the bi-abundance model B at a viewing angel of 8 = 60 deg. 9 = deg is for face on. The units of the color bars are 
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Figure 15. Comparison of the predicted and observed line fluxes for ORLs (left panel) and CELs (right panel) in the bi-abundance 
model B. Black and gray parts represent contributions from the normal and the cold H-deficient components, respectively. 



is not decreased compared to prediction. We also 
want to point out that in the case of the escaping 
He I Lyman photons, the discrepancy factors for the 
np 1 P-2s 1 S series (n = 3,4,5) will decrease with in- 
creasing n values, which is also not seen in the Orion 
Nebula. In their non-Case-B, self-consistent "model 
III", Porter et al. (2007) attributed the discrepan- 
cies to "inaccurate reddening corrections". Another 
plausible mechanism is that He I Lyman photons 
are absorbed by hydrogen atoms or dust grains, as 
proposed by Liu et al. (2001). We conclude that while 
CaseB recombination may still remain a good approxima- 
tion for the He I singlet lines in some PNe (e.g. NGC 7027; 
Zhang et al. 2005b), significant departures from this approx- 
imation may occur in others, and the effects can be quite 
significant for the 2s 1 S-np 1 P and 2p 1 P-ns 1 S series. 

We deduce the extinction coefficient c(H/3) by com- 
paring the observed Hi and He II line ratios to theoret- 
ical values, assuming a constant electron temperature and 
density. In the presence of a cold, H-deficient component, 



such as in the bi-abundance model, the extinction coefficient 
thus derived would be overestimated, leading to the fluxes 
of lines of short wavelengths being systematically overesti- 
mated while those of long wavelengths underestimated. This 
effect is nevertheless small, less than 5 per cent in the current 
bi-abundance model. 

4.2 T fluctuations versus the bi-abundance 
model 

T c fluctuations have been proposed to explain the 
CEL/ORL abundance discrepancies as well as the differ- 
ences between the electron temperatures derived from the 
[O ill] CELs and from the hydrogen Balmer discontinuity 
(e.g. Peimbert & Peimbert, 2006). However, chemically ho- 
mogeneous models predict a very small t 2 . Using the for- 
mula of Peimbert (1967) and the values of T ([0 ill]) and 
T C ([H i]) measured by L2000, we find t 2 = 0.065, much larger 
than the values predicted by chemically homogeneous mod- 
els (t 2 ^ 0.008; Table 2). The bi-abundance model yields a 
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larger value (t 2 ~ 0.014), but still significantly lower than 
the 0.065 derived from the measured T e ([0 ill]) and T e ([H i]). 
Note that the formalism of Peimbert (1967) is only applica- 
ble to situations of small amplitudes of temperature varia- 
tions. It will significantly overestimate t 2 in the case where 
the nebula is composed of two distinct components of very 
different temperatures (Zhang et al. 2007). By satisfactorily 
reproducing observations of the [O ill] CELs and the hydro- 
gen Balmer discontinuity, the bi-abundance model provides 
a natural explanation for the discrepancy between T c ([0 ill]) 
andT e ([Hi]). 

4.3 Properties of the cold component 

Table 4 summarizes the properties of the H-deficient knots 
derived from the bi-abundance model. The model shows that 
the knots are fully ionized and optically thin, with T a of 
800 K, a hydrogen number density of 4,000 cm~ 3 , a He/H 
abundance ratio of 0.5, and CNONe abundances relative to 
hydrogen about 40-100 times higher than the "normal" com- 
ponent. The mass of this cold H-deficient component is only 
0.0031 Mq (~ 3 Jupiter masses), about 1 per cent of the 
total mass of the nebula. 

With much enhanced metal abundances compared to 
the "normal" component, the cooling of the cold compo- 
nent is dominated by the IR fine-structure lines from C, N, 
O and Ne ions. These fine-structure lines have different ex- 
citation temperatures, typically much lower than those of 
optical/UV CELs. Their emissivities are thus less sensitive 
to temperature than optical/UV CELs. These IR lines also 
have relatively low, yet varying critical densities. Therefore, 
depending on the local physical conditions, their contribu- 
tions to the cooling vary. It follows that the electron den- 
sity plays an important role in setting the equilibrium tem- 
peratures of the cold component by determining which IR 
lines are dominant coolants. At low electron densities, the 
[O ill] 52 and 88 /jm lines are the dominant coolants, whereas 
at higher densities, the [Ne n] 12.8 pm and [Ne ill] 15.6, 
36 /J,m lines, owing to their higher critical densities and 
higher excitation temperatures than the [O ill] lines, start 
to dominate the cooling. Because the [Ne n] and [Ne ill] IR 
lines have higher excitation temperatures than the [O ill] IR 
lines, the equilibrium temperatures in high-density knots are 
higher than those in low-density ones. For example, as hy- 
drogen number density decreases from 10,000 cm~ 3 to 1,000 
cm~ 3 , the equilibrium temperature drops from 1,070 K to 
460 K. It follows that both the density and temperature 
contrasts between the cold component and the normal com- 
ponent are expected to grow as the nebula expands, assum- 
ing that the two components are in pressure equilibrium. 
As a consequence, one expects the adf value to increase as 
the nebula evolves, consistent with what is observed (e.g. 
Zhang et al. 2004). The equilibrium temperatures in the 
H-deficient knots are also affected by the He and heavy el- 
ement abundances in that the photoionization of He domi- 
nates the heating whereas the CNONe IR lines control the 
cooling. We find that an increase of 60 per cent in He abun- 
dance or a decrease of 20 per cent in CNONe abundances 
leads to an increase of 200 K in the equilibrium tempera- 
ture. It should be mentioned that, because of the lack of 
suitable diagnostic tools, densities of the H-deficient inclu- 
sions are not well constrained in the current models. This 
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Figure 16. Pressure ratios of cold H-dcficicnt knots and their 
surrounding gas in the bi-abundance model. 



introduces some uncertainties in the estimate of the total 
mass of the cold component. It is also worth pointing out 
that in the bi-abundance model, electron temperatures de- 
duced from ratios of IR CELs to optical CELs, such as the 
[Ne ill] 15.6/im/A3869 ratio, are expected to be lower than 
the value derived from the [O ill] optical line ratio, consis- 
tent with the observations (e.g. Bernard-Salas et al. 2004). 

We compare the pressures of the cold inclusions and 
their surrounding gas in Fig. 16. On average, the two compo- 
nents have a pressure ratio P(cold) /P(normal) = 0.53±0.77, 
which means that the H-deficient knots have a large range of 
pressure relative to the surrounding gas. This departure 
from pressure equilibrium could be a consequence of 
some simplifying assumptions of the model. In the 
model we assume that the H-deficient knots have 
the same hydrogen number density, and the main 
nebula has a smoothly varying density distribution 
except the equatorial torus. However, the departure 
is not deleterious to our conclusions. 

The H-deficient knots are thermally stable. If the elec- 
tron temperature increases, the enhanced thermal pressure 
will force the knots to expand and thus decrease the density, 
leading to enhanced cooling that will reduce the tempera- 
ture. This property may help survival of H-deficient knots 
during the PN evolution phase. 

In the current bi-abundance model, the cold H-deficient 
knots are distributed in the nebular inner region (within 10 
arcsec) in order to reproduce the ORL surface brightness dis- 
tributions yielded by long-slit observations (L2000) . It seems 
to be a common characteristic of PNe that ORL strengths 
peak at the nebular centre, e.g., PN NGC 7009 (Luo et al. 
2001; Tsamis et al. 2008), PN NGC 6720 (Garnett & Din- 
erstein 2001). Our modelling also shows that in order to 
match the strength of the [Ne n] 12.8 /im line, which orig- 
inates mainly from the cold component, the knots must be 
located relatively close to the central star in order to avoid 
over-production of Ne + ions. It is unclear whether this con- 
centration of H-deficient inclusions near the nebular centre 
is a consequence of their lower expansion velocity compared 
to the main diffuse gas or whether they form from a later 
evolutionary phase of the central star. 

Given their high densities, high CNONe abundances 
and low ionization degrees, the cold knots absorb hard ion- 
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Table 4. Properties of the H-dcficicnt inclusions in model B. 



(To) (K) 


815 


(JV H > (cm- 3 ) 


4000 


(N e ) (cm" 3 ) 


6 680 


M (M Q ) 


0.0031 


Mass fraction 


1.3 % 


Number of cells 


872 


Cell size 


0'.'25 x 0"25 x 0'.'33 


Filling factor 


0.002 


He/H 


0.50 


C/H 


0.01177 


N/H 


0.0150 


O/H 


0.0440 


Ne/H 


0.0177 



izing photons more efficiently than the diffuse gas, and 
thus weaken and soften the stellar radiation fields passing 
through them. As a result, the radial shadow tails behind 
the H-deficient knots have lower temperatures and ioniza- 
tion degrees than their nearby "unshielded" gas. These ef- 
fects are visible in Figs. 12 and 13. The presence of these 
shadow regions is however unlikely to have major effects on 
ORLs and low-ionization CELs or on t 2 values. 

Finally, we point out that the current investigation of 
the properties of H-deficient inclusions can be much im- 
proved with better IR observations. Restricted by the com- 
puter memory, the knots are not spatially resolved at all 
in the current modelling. This rules out the possibility to 
study structure and properties of optically thick H-deficient 
knots as well as effects of knots smaller than the model res- 
olution limit. In addition, the current work does not take 
into account the effects of dust grains, a potential heating 
source (e.g. Stasihska & Szczerba 2001). Given the differ- 
ent physical and chemical conditions, the dust grains in H- 
deficient knots presumably differ in properties from those in 
the "normal" gas. Future modelling that takes into account 
the effects of dust grains on nebular ionization and thermal 
structures may be worthwhile. 

4.4 Size of H-deficient knots 

The grainy H/3 image predicted by model B (Fig. 14) is 
not consistent with the HST observation, suggesting that 
size of the postulated H-deficient inclusions must be smaller 
than adopted in this model. In order to investigate possible 
effects of assumed knot size on model results, we have con- 
structed another bi-abundance model of higher resolution 
using a grid of 96 3 cells (model B'). For this new model, all 
the input parameters are the same as those used in model 
B except that the physical size of knots (cells) is reduced by 
half and the number of knots is increased by a factor of 8. 
Fig. 17 shows the electron density and temperature struc- 
tures predicted by model B'. The resultant ORL and CEL 
images are displayed in Fig. 18. The predicted total line 
fluxes as well as the contributions from the two components 
are presented in the last three columns of Table 1. For most 
lines, the differences between models B and B' are small. 
The cold H-deficient knots in model B' have slightly higher 
ionization degrees than those in model B due to their smaller 
physical size in model B'. There are however significant dif- 
ferences between the two models for the ionic fractions of 



high-ionization ions such as C 3+ , N 3+ , O s+ and Hc 2+ as 
well as low ionization ions such as Ne + . As a result, the 
predicted fluxes of lines from those ions vary dramatically. 
For example, the [Ne n] 12.8 /im line flux emitted by the 
cold H-deficient knots in model B' decreases by 50 per cent 
compared to model B, while the fluxes of C ill, N ill ORLs 
emitted by the cold H-deficient knots increase by 50 per cent, 
and the fluxes of He II lines and the [O iv] 25.9 /im line are 
almost doubled. These results imply that lines such as the 
[Ne n] 12.8 /itm CEL and C ill, N ill ORLs can potentially be 
used to constrain the size of the cold H-deficient inclusions. 
This is however complicated by the fact that intensities of 
these lines are also affected by the electron density in the 
cells as well as their spatial distribution. 

In an attempt to better constrain the size of the cold 
knots, the STIS spectra obtained by us have been used to 
deduce the surface brightness distribution of the C II A4267 
line, the strongest ORL detected by STIS. The spectra are 
unfortunately quite noisy and suffer from severe cosmic ray 
contamination. Nevertheless, we find the observed surface 
brightness distribution across the nebula is smoother than 
predicted by model B', suggesting that the postulated cold 
H-deficient inclusions must have sizes smaller than 1/6" (i.e. 
250 AU at a distance of 1.5 kpc). 

4.5 Elemental abundances 

In Table 5, we list the chemical abundances of NGC 6153 
derived from the bi-abundance model B and those in the 
literature obtained from traditional empirical analyses. For 
comparison, we also give the solar abundances presented by 
Asplund et al. (2009). Table 5 shows that the heavy element 
abundances deduced from optical/UV CELs, which are sys- 
tematically lower than values yielded by infrared CELs (Pot- 
tasch et al. 2003), are in excellent agreement with those of 
the "normal" component; whereas the heavy element abun- 
dances deduced from infrared CELs are close to the aver- 
age values of the whole nebula, and are much lower than 
the abundances in the cold component. The helium abun- 
dance published in the literature, deduced from the empiri- 
cal method assuming a homogeneous composition, has been 
significantly overestimated due to the contamination of cold 
knots (see also Pequignot et al. 2002). If this is also the case 
for H II regions, the primordial helium abundance deduced 
from the standard analyses of helium recombination lines of 
metal-poor galaxies will be overestimated. The He, C, N, O 
and Ne abundances in the cold knots are respectively 5, 55, 
40, 80 and 100 times higher than those in the main nebula. 

4.6 Possible origins of the postulated H-deficient 
inclusions 

Given that adf values measured for PNe are always larger 
than unity, H-deficient inclusions are presumably a gen- 
uine feature of PNe. Their presence in PNe is however not 
predicted by current theories of stellar evolution. While it 
has been proposed (Iben et al. 1983) that an evolved star 
undergoing a very late helium flash (the so called "born- 
again" PNe) may harbor H-deficient material, the central 
stars of most PNe exhibiting large adf 's are found not to be 
H-deficient. Hydrogen-deficient knots have been clearly de- 
tected in the two "born-again" PNe, Abell30 and Abell58. 
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Three-dimensional structures of Af c (left panel) and T e (right panel) of model B' that used a higher spatial resolution of 96 3 
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Figure 18. Projected monochromatic images of H/3 (left panel), [O III] A5007 (central panel) and O II A4649 (right panel) of model B' 



at a viewing angel of 
respectively. 



60 deg. The units of the color bar are 10 13 , 10 12 and 10 15 ergs cm 2 s 1 arcsec 2 for the three images, 



Table 5. Elemental abundances of NGC 6153. 



Elem. 


L2000 a 


P03 6 


B n 


B c 


B 


Solar c 


He 


0.14 


0.14 


0.10 


0.50 


0.102 


0.085 


C(-4) 


2.8 


6.8 


3.2 


177 


3.88 


2.692 


N(-4) 


2.3 


4.8 


3.8 


150 


4.37 


0.676 


0(-4) 


5.0 


8.3 


5.63 


440 


7.33 


4.898 


Ne(-4) 


1.7 


3.1 


1.76 


177 


2.44 


0.851 


Mg(-5) 






3.8 


12.1 


3.83 


3.981 


Si(-5) 






3.5 


11.3 


3.53 


3.236 


S(-5) 


1.6 


1.9 


1.75 


5.16 


1.76 


1.318 


Cl(-7) 


4.2 


5.6 


2.35 


10.1 


2.38 


3.162 


Ar(-6) 


2.7 


8.5 


2.9 


11.5 


2.93 


2.512 



a from L2000 for optical/UV CELs. 
b from Pottasch et al. (2003). 
c from Asplund et al. (2009). 



On the other hand, Wesson et al. (2003, 2008) found that the 
H-deficient knots in Abell30 and Abell 58 are oxygen-rich, 
in contrast to the expectation of the "born-again" scenario. 
Their chemical composition has more in common with neon 
novae than with Sakurai's Object. The latter is believed to 
have recently experienced a final He-shell flash. Lau et al. 
(in preparation) explore the possibility of binary evolution to 
account for the observed high oxygen and neon abundances 
of the H-deficient knots in Abell 58, assuming that the knots 
are from the ejecta of a neon nova explosion. They consider a 
number of scenarios but none fit perfectly with the observa- 
tions. Their best scenario consists of binary stars where the 
primary ONe white dwarf provided the neon nova explosion 
and the secondary asymptotic giant branch star evolved into 
the PN. This scenario requires very fine-tuned initial condi- 
tions, particularly the initial separation, to make the neon 
nova occur just after the final flash of the secondary asymp- 
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totic giant branch companion. Therefore, this scenario does 
not predict a large number of PNe with H-deficient knots. 

Another conjecture is that the H-deficient inclusions are 
evaporating solid bodies, such as metal-rich planetesimals 
formed in planetary disks (Liu et al. 2006). This scenario 
has been subsequently investigated by Henney & Stasihska 
(2010), who show that the destruction of solid bodies during 
the short PN phase itself is not able to provide enough metal- 
rich material to explain the observed CEL/ORL abundance 
discrepancies, due to a low sputtering rate. On the other 
hand, the sublimation of solid bodies during the final stages 
of the asymptotic giant branch phase, which is much longer 
than the PN phase, might provide enough metal-rich ma- 
terial and serve as a possible mechanism causing the abun- 
dance discrepancies. 

Although other explanations to the (moderate) ADF 
of HIIRs are still possible, the H-deficient inclusions may 
also be present in H II regions (Tsamis et al. 2003; Tsamis 
& Pequignot, 2005). They could be supernovae ejecta in the 
form of metal-rich droplets that fall back into the interstellar 
medium after a long journey lasting about 10 8 yr (Stasihska 
et al. 2007). They might also originate from evaporating 
proto-planetary disks around newly formed stars, such as 
the proplyds discovered in the Orion Nebula. 



5 SUMMARY 

We have constructed detailed three-dimensional photoion- 
ization models of NGC 6153 of spherical or ellipsoidal ge- 
ometry, with and without H-deficient inclusions. The model 
results are compared with spectroscopic observations from 
the UV, optical to the far IR as well as narrow-band imaging 
observations. In those models, the main nebula is assumed to 
be chemically homogeneous and of a "normal" composition. 
The spherical model can reproduce the strengths of hydro- 
gen and helium recombination lines and of high excitation 
CELs. In order to explain the strengths of low-excitation 
CELs, we have considered an ellipsoidal shell where the den- 
sity decreases from the equator to the poles, supplemented 
with an equatorial torus having the same chemical compo- 
sition but of a higher density. The chemically homogeneous 
ellipsoidal model with a torus reproduces the strengths of 
all CELs within the measurement uncertainties, but fails 
by an order of magnitude for heavy element ORLs. In all 
cases, chemically homogeneous models yield small temper- 
ature fluctuations, and are therefore incapable to explain 
the large discrepancy between T c ([0 III]) and T (H I). In 
contrast, a bi-abundance model, incorporating small-scale 
chemical inhomogencities in the form of a small amount of 
metal-rich inclusions, not only reproduces the strengths of 
heavy element ORLs and the hydrogen Balmer discontinu- 
ity, but also improves the overall fitting of CELs. 

In our best-fit bi-abundance model, the CNONe abun- 
dances in the metal-rich knots are enhanced by 40-100 times 
compared to those of the main nebula. The presence of those 
metal-rich knots increases the average metal content of the 
whole nebula by 30 per cent. Optical/UV CELs arise mainly 
from the main diffuse nebula of "normal" composition and 
thus probe the physical conditions and chemical abundances 
of this component. IR CELs originate from both compo- 
nents, and thus yield elemental abundances higher than op- 



tical/UV CELs but lower than ORLs. In addition, temper- 
atures deduced from IR CELs tend to be lower than values 
deduced from optical/UV CELs. For He I lines, although 
a large fraction of their fluxes arises from the H-deficient 
component, given the relatively low abundance enhancement 
factor of about 5, the helium abundance of the normal com- 
ponent remains a fair representation of the average value 
of the whole nebula. Nevertheless, in the presence of H- 
deficient inclusions, the traditional analysis assuming a sin- 
gle uniform chemical composition for the whole nebula, will 
greatly overestimate the helium abundances of the normal 
component. One should bear in mind that in the current 
work we were unable to separate helium emission originat- 
ing from the two components to obtain reliable estimates 
of the helium abundance in each of the two components. 
Indeed, we have assumed He/H = 0.1 for the normal com- 
ponent. The collisional excitation dominant He I A 10830 line 
can, in principle, be used to estimate the helium abundance 
in the normal component. Its usage is however hampered 
by observational and theoretical uncertainties. The electron 
temperature inside the metal-rich knots is as low as 800 K. 
The main coolants are IR fine-structure lines from heavy el- 
ement ions, including the [O in] 52, 88 /im, the [Ne ill] 15.6, 
36 fim and [Ne n] 12.8 pm lines. The metal-rich knots have 
a high density, and are roughly in pressure equilibrium with 
the surrounding "normal" gas. The metal-rich knots weaken 
and soften the UV radiation field. The cometary tails of 
"normal" gas behind the knots have relatively low temper- 
atures and ionization degrees. 

We show that the presence of metal-rich knots may lead 
to significant errors in nebular properties derived from the 
traditional methods assuming a single nebular component. 
We discuss physical and chemical properties of the metal- 
rich knots in NGC 6153 as well as their possible origins. To 
better constrain the nature of metal-rich knots in PNe, new 
observations and techniques are useful. For example, both 
imaging and spectroscopic observations of IR fine-structure 
lines at high spatial resolution and sensitivity will certainly 
shed new insight into the properties of these metal-rich 
knots. In future work, we plan to develop more comprehen- 
sive methods to probe these cold knots by involving more 
observational features (e.g. He I discontinuities; Zhang et 
al. 2005a, 2009) and up-to-date recombination coefficients 
of ORLs for low-temperature plasma (Fang et al. in prepa- 
ration). We will also construct detailed three-dimensional 
models for more PNe exhibiting large adf's with resolved 
metal-rich knots, such as M 1-42. These are the subjects of 
further papers. 
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